% This function returns EV for an individual who chooses the default rate
% Note: m0 is the amount by which income is augmented by the equalizing compensation (currently zero each time function is called), which can affect the contribution rate choice. In our final specifications, we keep the contribution rate fixed as we vary equalizing compensation (see footnote 24).

function[value]=ev_optin_int(m0,alpha,rho,t,m,maxmatch,d,maxcont,mm)

% Define virtual income above the kink
Zbar=1+(1-t).*maxmatch.*(1+m).*(1-(1/(1+m)));
d_matched = d + m.*min(d,maxmatch);

% Compute EV. There are four cases to consider here, because x*(theta)
% and the default can both be above of below the maximum match rate.
ev_ifxdbelowM = (((d_matched+alpha)./(xstar_ev(m0,alpha,rho,t,m,Zbar,maxmatch,maxcont)+alpha)).^rho).*...
    (1-((1-t)./(1+m)).*d_matched)+(((1-t)./(1+m)).*xstar_ev(m0,alpha,rho,t,m,Zbar,maxmatch,maxcont))- 1 +mm.*d_matched ;

ev_ifxdaboveM = (((d_matched+alpha)./(xstar_ev(m0,alpha,rho,t,m,Zbar,maxmatch,maxcont)+alpha)).^rho).*...
    (Zbar-((1-t)).*d_matched)+(((1-t)).*xstar_ev(m0,alpha,rho,t,m,Zbar,maxmatch,maxcont))- Zbar +mm.*d_matched ;

ev_ifxbelowdaboveM = (((d_matched+alpha)./(xstar_ev(m0,alpha,rho,t,m,Zbar,maxmatch,maxcont)+alpha)).^rho).*...
    (Zbar-((1-t)).*d_matched)+(((1-t)./(1+m)).*xstar_ev(m0,alpha,rho,t,m,Zbar,maxmatch,maxcont))- 1 +mm.*d_matched ;

ev_ifxabovedbelowM = (((d_matched+alpha)./(xstar_ev(m0,alpha,rho,t,m,Zbar,maxmatch,maxcont)+alpha)).^rho).*...
    (1-((1-t)./(1+m)).*d_matched)+(((1-t)).*xstar_ev(m0,alpha,rho,t,m,Zbar,maxmatch,maxcont))- Zbar + mm.*d_matched ;

value = ev_ifxdbelowM.*(xstar_ev(m0,alpha,rho,t,m,Zbar,maxmatch,maxcont)<=maxmatch.*(1+m)).*(d<=maxmatch)...
    + ev_ifxdaboveM.*(xstar_ev(m0,alpha,rho,t,m,Zbar,maxmatch,maxcont)>maxmatch.*(1+m)).*(d>maxmatch)...
    + ev_ifxbelowdaboveM.*(xstar_ev(m0,alpha,rho,t,m,Zbar,maxmatch,maxcont)<=maxmatch.*(1+m)).*(d>maxmatch)...
    + ev_ifxabovedbelowM.*(xstar_ev(m0,alpha,rho,t,m,Zbar,maxmatch,maxcont)>maxmatch.*(1+m)).*(d<=maxmatch);




    % x*(theta,m0) function
    function[value]=xstar_ev(m1,alpha,rho,t,m,Zbar,maxmatch,maxcont)
    xstar_ifbelowM = (rho/(1+rho)).*((1+m).*(1+m1)./(1-t))-alpha./(1+rho);
    xstar_ifaboveM = (rho/(1+rho)).*((Zbar+m1)./(1-t))-alpha./(1+rho);
    xstar = xstar_ifbelowM.*(xstar_ifbelowM<=maxmatch.*(1+m)) + xstar_ifaboveM.*(xstar_ifaboveM>maxmatch.*(1+m))...
        + maxmatch.*(1+m).*(xstar_ifbelowM>maxmatch.*(1+m)).*(xstar_ifaboveM<=maxmatch.*(1+m));
    xstar=max(xstar,0);
    
    % Now, we account for the fact that the personal contribution,
    % x_personal, must be a discrete value (i.e. 1%, 2%, etc.)
    x_personal_opt=min(maxmatch.*(1+m),xstar)./(1+m)+max(0,xstar-(maxmatch.*(1+m)));
    
    % Calculate discrete % contributions above and below optimal personal
    % contribution. (Can't contribute < 0)
    x_personal_top=ceil(x_personal_opt.*100)./100;
    x_personal_top=max(x_personal_top,0);
    
    x_personal_bot=floor(x_personal_opt.*100)./100;
    x_personal_bot=max(x_personal_bot,0);
    
    % Can't contribute > maxcont
    x_personal_top=min(maxcont,x_personal_top);
    x_personal_bot=min(maxcont,x_personal_bot);
    
    % So the choice, x, would actually be one of these two discrete values
    x_top=x_personal_top+(m.*min(x_personal_top,maxmatch));
    x_bot=x_personal_bot+(m.*min(x_personal_bot,maxmatch));
    
    % The employee chooses the round number that gives the highest utility
    utility_top = ((rho.*log(x_top+alpha) + log(1+m1-((1-t)./(1+m)).*x_top)).*...
        (x_top<=maxmatch.*(1+m))... 
        +((rho.*log(x_top+alpha)+log(Zbar+m1-(1-t).*x_top)).*((x_top>maxmatch.*(1+m)))));
    utility_bot = ((rho.*log(x_bot+alpha) + log(1+m1-((1-t)./(1+m)).*x_bot)).*...
        (x_bot<=maxmatch.*(1+m)))... 
        +((rho.*log(x_bot+alpha)+log(Zbar+m1-(1-t).*x_bot)).*((x_bot>maxmatch.*(1+m))));
    
    value=(xstar>=0).*( (x_top).*(utility_top>=utility_bot) + (x_bot).*(utility_top<utility_bot));
    end

end